###Do Partisan Types Stop at the Water's Edge? 
#Partisan Types 2.R
#Josh Kertzer
#Last modified: February 13, 2020

### Note: this file runs the analysis in the supplementary appendix - run Partisan Types 0.R and Partisan Types 1.R first
### All of the following analyses were carried out using R version 3.5.2 on a 3 GHz Intel Xeon W iMac Pro running macOS Mojave 10.14.6, and replicated using R version 3.5.2 on a 2 GHz Intel Core i7 Macbook Pro running macOS Mojave 10.14.5.

#### Appendix 1.3: Sampling methodology

## Table 5: Sample characteristics

tab2014 <- c(round(prop.table(table(dat1$male)), digits=3), round(prop.table(table(dat1$age_1)), digits=3)[2], round(prop.table(table(dat1$age_2)), digits=3)[2], round(prop.table(table(dat1$age_3)), digits=3)[2], round(prop.table(table(dat1$age_4)), digits=3)[2], round(prop.table(table(dat1$educ_1)), digits=3)[2], round(prop.table(table(dat1$educ_2)), digits=3)[2], round(prop.table(table(dat1$educ_3)), digits=3)[2], round(prop.table(table(dat1$educ_4)), digits=3)[2], round(prop.table(table(dat1$white)), digits=3)[2], round(prop.table(table(dat1$pid3)), digits=3))

tab2018 <- c(round(prop.table(table(dat2$male)), digits=3), round(prop.table(table(dat2$age_1)), digits=3)[2], round(prop.table(table(dat2$age_2)), digits=3)[2], round(prop.table(table(dat2$age_3)), digits=3)[2], round(prop.table(table(dat2$age_4)), digits=3)[2], round(prop.table(table(dat2$educ_1)), digits=3)[2], round(prop.table(table(dat2$educ_2)), digits=3)[2], round(prop.table(table(dat2$educ_3)), digits=3)[2], round(prop.table(table(dat2$educ_4)), digits=3)[2], round(prop.table(table(dat2$white)), digits=3)[2], round(prop.table(table(dat2$pid3)), digits=3))

dat2$Dem <- recode(dat2$pid3, binarize=TRUE, x=1)
dat2$Ind <- recode(dat2$pid3, binarize=TRUE, x=2)
dat2$Rep <- recode(dat2$pid3, binarize=TRUE, x=3)
dat2$female <- recode(dat2$male, binarize=TRUE, x=0)

#Remove observations without weights
dat2a <- dat2[which(!is.na(dat2$weight)),]
tab2018w <- c(round(weighted.mean(dat2a$female, dat2a$weight, na.rm=TRUE), digits=3), round(weighted.mean(dat2a$male, dat2a$weight, na.rm=TRUE), digits=3), round(weighted.mean(dat2a$age_1, dat2a$weight, na.rm=TRUE), digits=3), round(weighted.mean(dat2a$age_2, dat2a$weight, na.rm=TRUE), digits=3), round(weighted.mean(dat2a$age_3, dat2a$weight, na.rm=TRUE), digits=3), round(weighted.mean(dat2a$age_4, dat2a$weight, na.rm=TRUE), digits=3), round(weighted.mean(dat2a$educ_1, dat2a$weight, na.rm=TRUE), digits=3), round(weighted.mean(dat2a$educ_2, dat2a$weight, na.rm=TRUE), digits=3), round(weighted.mean(dat2a$educ_3, dat2a$weight, na.rm=TRUE), digits=3), round(weighted.mean(dat2a$educ_4, dat2a$weight, na.rm=TRUE), digits=3), round(weighted.mean(dat2a$white, dat2a$weight, na.rm=TRUE), digits=3), round(weighted.mean(dat2a$Dem, dat2a$weight, na.rm=TRUE), digits=3), round(weighted.mean(dat2a$Ind, dat2a$weight, na.rm=TRUE), digits=3), round(weighted.mean(dat2a$Rep, dat2a$weight, na.rm=TRUE), digits=3))

names(tab2014) <- names(tab2018) <- names(tab2018w) <- c("Female", "Male", "Age: 18-24", "Age: 25-44", "Age: 45-64", "Age: 65+", "Education: High school or less", "Education: Some college", "Education: College/university", "Education: Graduate/professional school", "White", "Democratic", "Independent", "Republican")

stargazer(cbind(tab2014, tab2018, tab2018w))

#### Appendix 2.1: Raw distributions of second-order beliefs

### Figure 1: Raw distributions of responses: study 1
## Study 1 ##
#First, use s.content2 to combine both/neither into middle category
content1 <- rescale(as.data.frame(apply(stereo1,2,s.content2, itemLevel=TRUE)))

#Now plot distributions of partisan types
vector.m <- na.omit(melt(content1))
levels(vector.m$variable) <- c("Stop piracy", "Stop piracy (uni)", "Stop piracy (multi)", "Enviro sanction", "Enviro sanction (uni)", "Enviro sanction (multi)", "Troops oil", "Troops oil (uni)", "Troops oil (multi)", "Troops humanitarian", "Troops humanitarian (uni)", "Troops humanitarian (multi)", "Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion")

vector.m$variable  <- relevel(vector.m$variable, "Stop piracy (multi)")
vector.m$variable  <- relevel(vector.m$variable, "Stop piracy (uni)")
vector.m$variable  <- relevel(vector.m$variable, "Stop piracy")
vector.m$variable  <- relevel(vector.m$variable, "Enviro sanction (multi)")
vector.m$variable  <- relevel(vector.m$variable, "Enviro sanction (uni)")
vector.m$variable  <- relevel(vector.m$variable, "Enviro sanction")
vector.m$variable  <- relevel(vector.m$variable, "Troops oil (multi)")
vector.m$variable  <- relevel(vector.m$variable, "Troops oil (uni)")
vector.m$variable  <- relevel(vector.m$variable, "Troops oil")
vector.m$variable  <- relevel(vector.m$variable, "Troops humanitarian (multi)")
vector.m$variable  <- relevel(vector.m$variable, "Troops humanitarian (uni)")
vector.m$variable  <- relevel(vector.m$variable, "Troops humanitarian")

vector.m$variable  <- relevel(vector.m$variable, "Decrease military spending")
vector.m$variable  <- relevel(vector.m$variable, "Increase military spending")
vector.m$variable  <- relevel(vector.m$variable, "Protectionist")
vector.m$variable  <- relevel(vector.m$variable, "Free trade")
vector.m$variable  <- relevel(vector.m$variable, "Isolationist")
vector.m$variable  <- relevel(vector.m$variable, "Interventionist")
vector.m$variable  <- relevel(vector.m$variable, "Anti-arms control")
vector.m$variable  <- relevel(vector.m$variable, "Pro-arms control")

vector.m$variable  <- relevel(vector.m$variable, "Pro-abortion")
vector.m$variable  <- relevel(vector.m$variable, "Anti-abortion")
vector.m$variable  <- relevel(vector.m$variable, "Pro-tax hike")
vector.m$variable  <- relevel(vector.m$variable, "Anti-tax hike")
vector.m$variable  <- relevel(vector.m$variable, "Pro-environment")
vector.m$variable  <- relevel(vector.m$variable, "Anti-environment")
vector.m$variable  <- relevel(vector.m$variable, "Pro-gun control")
vector.m$variable <- relevel(vector.m$variable, "Anti-gun control")

dev.new(height=8.5, width=12)
ggplot(vector.m, aes(x=value)) + geom_bar() + facet_wrap(~variable, ncol=7) + theme_bw() + labs(x="Distribution of responses by policy proposal", y="") + theme(axis.text.y=element_blank(), axis.ticks.y=element_blank()) + scale_x_continuous(breaks=c(0,0.5,1),labels=c("0"="Def.\nDem", "0.5"="Both/\nNeither","1"="Def.\nRep"))

### Figure 2: Raw distributions of responses: study 2

## Study 2 ##

#First, use s.content2 to combine both/neither into middle category
content2 <- rescale(as.data.frame(apply(stereo2,2,s.content2, itemLevel=TRUE)))

#Now plot distributions of partisan types
vector.m <- na.omit(melt(content2))
levels(vector.m$variable) <- c("Troops oil", "Troops oil (uni)", "Troops oil (multi)","Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion", "Anti-affirmative action", "Pro-affirmative action", "Anti-immigration", "Pro-immigration")

vector.m$variable  <- relevel(vector.m$variable, "Troops oil (multi)")
vector.m$variable  <- relevel(vector.m$variable, "Troops oil (uni)")
vector.m$variable  <- relevel(vector.m$variable, "Troops oil")

vector.m$variable  <- relevel(vector.m$variable, "Decrease military spending")
vector.m$variable  <- relevel(vector.m$variable, "Increase military spending")
vector.m$variable  <- relevel(vector.m$variable, "Protectionist")
vector.m$variable  <- relevel(vector.m$variable, "Free trade")
vector.m$variable  <- relevel(vector.m$variable, "Isolationist")
vector.m$variable  <- relevel(vector.m$variable, "Interventionist")
vector.m$variable  <- relevel(vector.m$variable, "Anti-arms control")
vector.m$variable  <- relevel(vector.m$variable, "Pro-arms control")
vector.m$variable  <- relevel(vector.m$variable, "Anti-immigration")
vector.m$variable  <- relevel(vector.m$variable, "Pro-immigration")

vector.m$variable  <- relevel(vector.m$variable, "Pro-affirmative action")
vector.m$variable  <- relevel(vector.m$variable, "Anti-affirmative action")
vector.m$variable  <- relevel(vector.m$variable, "Pro-abortion")
vector.m$variable  <- relevel(vector.m$variable, "Anti-abortion")
vector.m$variable  <- relevel(vector.m$variable, "Pro-tax hike")
vector.m$variable  <- relevel(vector.m$variable, "Anti-tax hike")
vector.m$variable  <- relevel(vector.m$variable, "Pro-environment")
vector.m$variable  <- relevel(vector.m$variable, "Anti-environment")
vector.m$variable  <- relevel(vector.m$variable, "Pro-gun control")
vector.m$variable <- relevel(vector.m$variable, "Anti-gun control")

dev.new(height=8.5, width=12)
ggplot(vector.m, aes(x=value)) + geom_bar() + facet_wrap(~variable, ncol=6) + theme_bw() + labs(x="Distribution of responses by policy proposal", y="") + theme(axis.text.y=element_blank(), axis.ticks.y=element_blank()) + scale_x_continuous(breaks=c(0,0.5,1),labels=c("0"="Def.\nDem", "0.5"="Both/\nNeither","1"="Def.\nRep"))

##### Appendix 2.2. Alternative measures of partisan types

### Figure 3: Stereotype prevalence, study 1

prevalence1 <- as.data.frame(ifelse(stereo1 <=4,1,ifelse(stereo1 >=5,0,NA)))
prevalence1$id <- seq_len(nrow(prevalence1))

B <- 1500
set.seed(43215)
prevalence.boot1 <- array(NA, dim=list(ncol(prevalence1)-1,B)) #Variable, bootstraps
for (i in seq_len(B)){
	for (j in 1:28){ #Variables
		resample <- sample(prevalence1$id, nrow(prevalence1), replace=T)
		temp <- prevalence1[resample,j]
		prevalence.boot1[j,i] <- mean(temp,na.rm=TRUE)
	}
}

prevalence.mat1 <- data.frame(V1=apply(prevalence.boot1,1,mean,na.rm=TRUE), V2=apply(prevalence.boot1,1,quantile,0.025,na.rm=TRUE), V3=apply(prevalence.boot1,1,quantile,0.975,na.rm=TRUE), V4=c(1:28))
prevalence.mat1$V5 <- "Domestic"
prevalence.mat1$V5[c(1:12,13:14,17:20,21:22)] <- "Foreign" 
prevalence.mat1$V6 <- factor(c("Stop piracy", "Stop piracy (unilateral)", "Stop piracy (multilateral)", "Enviro sanction", "Enviro sanction (unilateral)", "Enviro sanction (multilateral)", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Troops humanitarian", "Troops humanitarian (unilateral)", "Troops humanitarian (multilateral)", "Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion"), levels=c("Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops humanitarian", "Troops humanitarian (unilateral)", "Troops humanitarian (multilateral)", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Enviro sanction", "Enviro sanction (unilateral)", "Enviro sanction (multilateral)", "Stop piracy", "Stop piracy (unilateral)", "Stop piracy (multilateral)"))
prevalence.mat1$V7 <- prevalence.mat1$V5
prevalence.mat1$V7[1:12] <- "Foreign (Approach)"
prevalence.mat1$V8 <- c(rep(1:4,each=3), rep(5:12,each=2))

dev.new(height=7, width=8.8)
ggplot(data=prevalence.mat1, aes(x=V1,y=V6)) + geom_point() + geom_segment(aes(x=V2, y=V6,xend=V3, yend=V6)) +theme_bw() + guides(color=FALSE, alpha=FALSE) + scale_color_grey() + facet_wrap(~V7, nrow=3, scales="free") +  labs(x="Stereotype prevalence", y="") + xlim(0,1)

### Figure 4: Stereotype prevalence, study 2

prevalence2 <- as.data.frame(ifelse(stereo2 <=4,1,ifelse(stereo2 >=5,0,NA)))
prevalence2$id <- seq_len(nrow(prevalence2))

B <- 1500
set.seed(43215)
prevalence.boot2 <- array(NA, dim=list(ncol(prevalence2)-1,B)) #Variable, bootstraps
for (i in seq_len(B)){
	for (j in 1:23){ #Variables
		resample <- sample(prevalence2$id, nrow(prevalence2), replace=T)
		temp <- prevalence2[resample,j]
		prevalence.boot2[j,i] <- mean(temp,na.rm=TRUE)
	}
}

prevalence.mat2 <- data.frame(V1=apply(prevalence.boot2,1,mean,na.rm=TRUE), V2=apply(prevalence.boot2,1,quantile,0.025,na.rm=TRUE), V3=apply(prevalence.boot2,1,quantile,0.975,na.rm=TRUE), V4=c(1:23))
prevalence.mat2$V5 <- "Domestic"
prevalence.mat2$V5[c(1:3,4:5,8:13,22:23)] <- "Foreign" 
prevalence.mat2$V6 <- factor(c("Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)","Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion", "Anti-affirmative action", "Pro-affirmative action", "Anti-immigration", "Pro-immigration"), levels=c("Pro-affirmative action", "Anti-affirmative action", "Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-immigration", "Anti-immigration", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)"))
prevalence.mat2$V7 <- prevalence.mat2$V5
prevalence.mat2$V7[1:3] <- "Foreign (Approach)"
prevalence.mat2$V8 <- c(rep(1,3), rep(2:11,each=2))

dev.new(height=7, width=8.8)
ggplot(data=prevalence.mat2, aes(x=V1,y=V6)) + geom_point() + geom_segment(aes(x=V2, y=V6,xend=V3, yend=V6)) +theme_bw() + guides(color=FALSE, alpha=FALSE) + scale_color_grey() + facet_wrap(~V7, nrow=3, scales="free") +  labs(x="Stereotype prevalence", y="")  + xlim(0,1)

### Figure 5: Stereotype content, alternate measure, study 1

#Alternate stereotype content measure (without both/neither)
s.content.alt <- function(x, itemLevel=TRUE){
	temp <- x
	#First, remove both/neither options
	temp[which(temp >=5)] <- NA
	if (itemLevel==TRUE){
		return(temp)
	}
	else if(itemLevel==FALSE){
		return(mean(temp, na.rm=TRUE))
	}
}

content1 <- rescale(as.data.frame(apply(stereo1,2,s.content.alt, itemLevel=TRUE)))
content1$id <- seq_len(nrow(content1))

B <- 1500
set.seed(43215)
content.boot1 <- array(NA, dim=list(ncol(content1)-1,B)) #Variable, bootstraps
for (i in seq_len(B)){
	for (j in 1:28){ #Variables
		resample <- sample(content1$id, nrow(content1), replace=T)
		temp <- content1[resample,j]
		content.boot1[j,i] <- mean(temp,na.rm=TRUE)
	}
}

content.mat1 <- data.frame(V1=apply(content.boot1,1,mean,na.rm=TRUE), V2=apply(content.boot1,1,quantile,0.025,na.rm=TRUE), V3=apply(content.boot1,1,quantile,0.975,na.rm=TRUE), V4=c(1:28))
content.mat1$V5 <- "Domestic"
content.mat1$V5[c(1:12,13:14,17:20,21:22)] <- "Foreign" 
content.mat1$V6 <- factor(c("Stop piracy", "Stop piracy (unilateral)", "Stop piracy (multilateral)", "Enviro sanction", "Enviro sanction (unilateral)", "Enviro sanction (multilateral)", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Troops humanitarian", "Troops humanitarian (unilateral)", "Troops humanitarian (multilateral)", "Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion"), levels=c("Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops humanitarian", "Troops humanitarian (unilateral)", "Troops humanitarian (multilateral)", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Enviro sanction", "Enviro sanction (unilateral)", "Enviro sanction (multilateral)", "Stop piracy", "Stop piracy (unilateral)", "Stop piracy (multilateral)"))
content.mat1$V7 <- content.mat1$V5
content.mat1$V7[1:12] <- "Foreign (Approach)"
content.mat1$V8 <- c(rep(1:4,each=3), rep(5:12,each=2))

dev.new(height=7, width=8.8)
ggplot(data=content.mat1, aes(x=V1,y=V6)) + geom_point() + geom_segment(aes(x=V2, y=V6,xend=V3, yend=V6)) +theme_bw() + guides(color=FALSE, alpha=FALSE) + scale_color_grey() + facet_wrap(~V7, nrow=3, scales="free") +  labs(x="Stereotype content", y="") + geom_line(aes(x=V1, y=V6, group=V8, alpha=0.8)) + xlim(0,1)
#Edit group ordering in the third panel manually

### Figure 6: Stereotype content, alternate measure, study 2

content2 <- rescale(as.data.frame(apply(stereo2,2,s.content.alt, itemLevel=TRUE)))
content2$id <- seq_len(nrow(content2))

B <- 1500
set.seed(43215)
content.boot2 <- array(NA, dim=list(ncol(content2)-1,B)) #Variable, bootstraps
for (i in seq_len(B)){
	for (j in 1:23){ #Variables
		resample <- sample(content2$id, nrow(content2), replace=T)
		temp <- content2[resample,j]
		content.boot2[j,i] <- mean(temp,na.rm=TRUE)
	}
}

content.mat2 <- data.frame(V1=apply(content.boot2,1,mean,na.rm=TRUE), V2=apply(content.boot2,1,quantile,0.025,na.rm=TRUE), V3=apply(content.boot2,1,quantile,0.975,na.rm=TRUE), V4=c(1:23))
content.mat2$V5 <- "Domestic"
content.mat2$V5[c(1:3,4:5,8:13,22:23)] <- "Foreign" 
content.mat2$V6 <- factor(c("Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)","Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion", "Anti-affirmative action", "Pro-affirmative action", "Anti-immigration", "Pro-immigration"), levels=c("Pro-affirmative action", "Anti-affirmative action", "Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-immigration", "Anti-immigration", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)"))
content.mat2$V7 <- content.mat2$V5
content.mat2$V7[1:3] <- "Foreign (Approach)"
content.mat2$V8 <- c(rep(1,3), rep(2:11,each=2))

dev.new(height=7, width=8.8)
ggplot(data=content.mat2, aes(x=V1,y=V6)) + geom_point() + geom_segment(aes(x=V2, y=V6,xend=V3, yend=V6)) +theme_bw() + guides(color=FALSE, alpha=FALSE) + scale_color_grey() + facet_wrap(~V7, nrow=3, scales="free") +  labs(x="Stereotype content", y="") + geom_line(aes(x=V1, y=V6, group=V8, alpha=0.8)) + xlim(0,1)
#Edit group ordering in the third panel manually

### Figure 7: Stereotype intensity, alternate measure, study 1

#Alternate stereotype intensity measure (without both/neither)
s.intensity.alt <- function(x, itemLevel=TRUE){
		temp <- x
		#First, remove both/neither the neutral midpoint
		temp[which(temp >=5)] <- NA 
		#Now, calculate intensity
		if(itemLevel==TRUE){
			return(abs(temp - 2.5)/1.5)
		}
		else if(itemLevel == FALSE){
			return(mean(abs(temp-2.5)/1.5, na.rm=TRUE))
		}
}

intensity1 <- rescale(as.data.frame(apply(stereo1,2,s.intensity.alt, itemLevel=TRUE)))
intensity1$id <- seq_len(nrow(intensity1))

B <- 1500
set.seed(43215)
intensity.boot1 <- array(NA, dim=list(ncol(intensity1)-1,B)) #Variable, bootstraps
for (i in seq_len(B)){
	for (j in 1:28){ #Variables
		resample <- sample(intensity1$id, nrow(intensity1), replace=T)
		temp <- intensity1[resample,j]
		intensity.boot1[j,i] <- mean(temp,na.rm=TRUE)
	}
}

intensity.mat1 <- data.frame(V1=apply(intensity.boot1,1,mean,na.rm=TRUE), V2=apply(intensity.boot1,1,quantile,0.025,na.rm=TRUE), V3=apply(intensity.boot1,1,quantile,0.975,na.rm=TRUE), V4=c(1:28))
intensity.mat1$V5 <- "Domestic"
intensity.mat1$V5[c(1:12,13:14,17:20,21:22)] <- "Foreign" 
intensity.mat1$V6 <- factor(c("Stop piracy", "Stop piracy (unilateral)", "Stop piracy (multilateral)", "Enviro sanction", "Enviro sanction (unilateral)", "Enviro sanction (multilateral)", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Troops humanitarian", "Troops humanitarian (unilateral)", "Troops humanitarian (multilateral)", "Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion"), levels=c("Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops humanitarian", "Troops humanitarian (unilateral)", "Troops humanitarian (multilateral)", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Enviro sanction", "Enviro sanction (unilateral)", "Enviro sanction (multilateral)", "Stop piracy", "Stop piracy (unilateral)", "Stop piracy (multilateral)"))
intensity.mat1$V7 <- intensity.mat1$V5
intensity.mat1$V7[1:12] <- "Foreign (Approach)"
intensity.mat1$V8 <- c(rep(1:4,each=3), rep(5:12,each=2))

dev.new(height=7, width=8.8)
ggplot(data=intensity.mat1, aes(x=V1,y=V6)) + geom_point() + geom_segment(aes(x=V2, y=V6,xend=V3, yend=V6)) +theme_bw() + guides(color=FALSE, alpha=FALSE) + scale_color_grey() + facet_wrap(~V7, nrow=3, scales="free") +  labs(x="Stereotype intensity", y="")  + xlim(0,1)

### Figure 8: Stereotype intensity, alternate measure, study 2

intensity2 <- rescale(as.data.frame(apply(stereo2,2,s.intensity.alt, itemLevel=TRUE)))
intensity2$id <- seq_len(nrow(intensity2))

B <- 1500
set.seed(43215)
intensity.boot2 <- array(NA, dim=list(ncol(intensity2)-1,B)) #Variable, bootstraps
for (i in seq_len(B)){
	for (j in 1:23){ #Variables
		resample <- sample(intensity2$id, nrow(intensity2), replace=T)
		temp <- intensity2[resample,j]
		intensity.boot2[j,i] <- mean(temp,na.rm=TRUE)
	}
}

intensity.mat2 <- data.frame(V1=apply(intensity.boot2,1,mean,na.rm=TRUE), V2=apply(intensity.boot2,1,quantile,0.025,na.rm=TRUE), V3=apply(intensity.boot2,1,quantile,0.975,na.rm=TRUE), V4=c(1:23))
intensity.mat2$V5 <- "Domestic"
intensity.mat2$V5[c(1:3,4:5,8:13,22:23)] <- "Foreign" 
intensity.mat2$V6 <- factor(c("Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)","Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion", "Anti-affirmative action", "Pro-affirmative action", "Anti-immigration", "Pro-immigration"), levels=c("Pro-affirmative action", "Anti-affirmative action", "Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-immigration", "Anti-immigration", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)"))
intensity.mat2$V7 <- intensity.mat2$V5
intensity.mat2$V7[1:3] <- "Foreign (Approach)"
intensity.mat2$V8 <- c(rep(1,3), rep(2:11,each=2))

dev.new(height=7, width=8.8)
ggplot(data=intensity.mat2, aes(x=V1,y=V6)) + geom_point() + geom_segment(aes(x=V2, y=V6,xend=V3, yend=V6)) +theme_bw() + guides(color=FALSE, alpha=FALSE) + scale_color_grey() + facet_wrap(~V7, nrow=3, scales="free") +  labs(x="Stereotype intensity", y="") +  xlim(0,1)

### Figure 9: Euclidean distance measure of stereotypicality (study 1)

#Generate reference vectors
uniform.dist <- c(1/6,1/6,1/3,1/6,1/6)
null.dist <- c(0,0,1,0,0)

content1 <- rescale(as.data.frame(apply(stereo1,2,s.content2, itemLevel=TRUE)))

#### Bootstrap Euclidean distance from null type
B <- 1500
set.seed(43215)
dist.1 <- array(NA, dim=list(ncol(content1),B,2))
for (i in seq_len(B)){
	resample <- sample(seq_len(nrow(content1)), nrow(content1), replace=T)
	temp <- content1[resample,]
	boot.val <- t(prop.table(apply(temp,2,table),2))
	dist.1[,i,1] <- as.matrix(dist(rbind(boot.val, uniform.dist)))[29,1:28]		
	dist.1[,i,2] <- as.matrix(dist(rbind(boot.val, null.dist)))[29,1:28]	
}

### Calculate quantities of interest

dist.boot1 <- data.frame(V1=c(apply(dist.1[,,1],1,mean,na.rm=TRUE), apply(dist.1[,,2],1,mean,na.rm=TRUE)), V2=c(apply(dist.1[,,1],1,quantile, 0.025,na.rm=TRUE), apply(dist.1[,,2],1,quantile, 0.025,na.rm=TRUE)), V3=c(apply(dist.1[,,1],1,quantile, 0.975,na.rm=TRUE), apply(dist.1[,,2],1,quantile, 0.975,na.rm=TRUE)), V4=c(1:28), V5=factor(rep(c("Uniform", "Null"),each=28), levels=c("Uniform", "Null")))

for.lab <- rep("Domestic",28)
for.lab[c(13:14,17:20,21:22)] <- "Foreign"
for.lab[c(1:12)] <- "Foreign (Approach)"
dist.boot1$V6 <- rep(for.lab,2)

dist.boot1$V4 <- factor(rep(c("Stop piracy", "Stop piracy (unilateral)", "Stop piracy (multilateral)", "Enviro sanction", "Enviro sanction (unilateral)", "Enviro sanction (multilateral)", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Troops humanitarian", "Troops humanitarian (unilateral)", "Troops humanitarian (multilateral)", "Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion"),2), levels=c("Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops humanitarian", "Troops humanitarian (unilateral)", "Troops humanitarian (multilateral)", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Enviro sanction", "Enviro sanction (unilateral)", "Enviro sanction (multilateral)", "Stop piracy", "Stop piracy (unilateral)", "Stop piracy (multilateral)"))

dev.new(height=7, width=8.8)
ggplot(data=dist.boot1, aes(x=V1,y=V4, col=V5)) + geom_point() + geom_segment(aes(x=V2, y=V4,xend=V3, yend=V4)) +theme_bw() + guides(alpha=FALSE, color=guide_legend(title="Reference\ndistribution")) + scale_color_grey() + facet_wrap(~V6, nrow=3, scales="free") +  labs(x="Euclidean distance from observed distribution to reference distribution", y="") + xlim(0,1.2) 

#### Figure 10: Euclidean distance measure of stereotypicality (study 2)

content2 <- rescale(as.data.frame(apply(stereo2,2,s.content2, itemLevel=TRUE)))

#### Bootstrap Euclidean distance from null type
B <- 1500
set.seed(43215)
dist.2 <- array(NA, dim=list(ncol(content2),B,2))
for (i in seq_len(B)){
	resample <- sample(seq_len(nrow(content2)), nrow(content2), replace=T)
	temp <- content2[resample,]
	boot.val <- t(prop.table(apply(temp,2,table),2))
	dist.2[,i,1] <- as.matrix(dist(rbind(boot.val, uniform.dist)))[24,1:23]		
	dist.2[,i,2] <- as.matrix(dist(rbind(boot.val, null.dist)))[24,1:23]	
}

### Calculate quantities of interest

dist.boot2 <- data.frame(V1=c(apply(dist.2[,,1],1,mean,na.rm=TRUE), apply(dist.2[,,2],1,mean,na.rm=TRUE)), V2=c(apply(dist.2[,,1],1,quantile, 0.025,na.rm=TRUE), apply(dist.2[,,2],1,quantile, 0.025,na.rm=TRUE)), V3=c(apply(dist.2[,,1],1,quantile, 0.975,na.rm=TRUE), apply(dist.2[,,2],1,quantile, 0.975,na.rm=TRUE)), V4=c(1:23), V5=factor(rep(c("Uniform", "Null"),each=23), levels=c("Uniform", "Null")))

for.lab <- rep("Domestic",23)
for.lab[c(1:3,4:5,8:13,22:23)] <- "Foreign"
for.lab[c(1:3)] <- "Foreign (Approach)"
dist.boot2$V6 <- rep(for.lab,2)

dist.boot2$V4 <- factor(rep(c("Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)","Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion", "Anti-affirmative action", "Pro-affirmative action", "Anti-immigration", "Pro-immigration"),2), levels=c("Pro-affirmative action", "Anti-affirmative action", "Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-immigration", "Anti-immigration", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)"))

dev.new(height=7, width=8.8)
ggplot(data=dist.boot2, aes(x=V1,y=V4, col=V5)) + geom_point() + geom_segment(aes(x=V2, y=V4,xend=V3, yend=V4)) +theme_bw() + guides(alpha=FALSE, color=guide_legend(title="Reference\ndistribution")) + scale_color_grey() + facet_wrap(~V6, nrow=3, scales="free") +  labs(x="Euclidean distance from observed distribution to reference distribution", y="") + xlim(0,1.2) 

#### Figure 11: Ternary plot of Euclidean distance measures

vector.1 <- t(prop.table(apply(content1, 2, table),2))
vector.2 <- t(prop.table(apply(content2, 2, table),2))

#Reference vectors
strong.dem <- c(1,0,0,0,0,0)
strong.rep <- c(0,0,0,1,0,0)
uniform.dist <- c(1/6,1/6,1/6,1/6,1/6,1/6)

dist.1d <- as.matrix(dist(rbind(vector.1, strong.dem)))[29,1:28] #Distance of each policy from a strong democratic reference point
dist.1r <- as.matrix(dist(rbind(vector.1, strong.rep)))[29,1:28] #Distance of each policy from a strong republican reference point
dist.1u <- as.matrix(dist(rbind(vector.1, uniform.dist)))[29,1:28] #Distance of each policy from uniform reference point
for.lab1 <- rep("Domestic",28)
for.lab1[c(1:12,13:14,17:20,21:22)] <- "Foreign" 

dist.2d <- as.matrix(dist(rbind(vector.2, strong.dem)))[24,1:23] #Distance of each policy from a strong democratic reference point
dist.2r <- as.matrix(dist(rbind(vector.2, strong.rep)))[24,1:23] #Distance of each policy from a strong republican reference point
dist.2u <- as.matrix(dist(rbind(vector.2, uniform.dist)))[24,1:23] #Distance of each policy from uniform reference point
for.lab2 <- rep("Domestic", 23)
for.lab2[c(1:3,4:5,8:13,22:23)] <- "Foreign"


tern.plot1 <- data.frame(Strong.Democrat=c(dist.1d,dist.2d), Strong.Republican=c(dist.1r,dist.2r), Uniform=c(dist.1u,dist.2u), for.lab=c(for.lab1,for.lab2))
dev.new(height=7, width=7)
ggtern(tern.plot1, aes(Strong.Democrat, Strong.Republican, Uniform, color=for.lab)) + geom_point(alpha=0.5) + theme_rgbw() + xlab("Distance:\nDemocratic") + ylab("Distance:\nRepublican") + zlab("Distance:\nUniform") + Larrowlab("Distance from pure\n Democratic type") + Tarrowlab("Distance from pure\n Republican type") + Rarrowlab("Distance from uniform\n distribution") + guides(color=guide_legend(title="Issue type"))
#Manually fix cropping

##### Figure 12: Stereotype content dropping "neither" category, study 1
s.content3 <- function(x, itemLevel=TRUE){
	temp <- x
	#First, re-order scale so that both is the neutral midpoint
	temp[which(temp==5)] <- 2.5 #Create new neutral midpoint for "both"
	temp[which(temp==6)] <- NA #Drop "neither"
	temp[which(temp>=3)] <- temp[which(temp>=3)]+1 #Incrase Republican responses by one
	temp[which(temp==2.5)] <- 3 #Now, move neutral midpoint up
	if (itemLevel==TRUE){
		return(temp)
	}
	else if(itemLevel==FALSE){
		return(mean(temp, na.rm=TRUE))
	}
}

s.contentN <- function(x, itemLevel=TRUE){
	temp <- x
	temp[which(temp<=5)] <- 0
	temp[which(temp==6)] <- 1
	if (itemLevel==TRUE){
		return(temp)
	}
	else if(itemLevel==FALSE){
		return(mean(temp, na.rm=TRUE))
	}	
}

content1 <- rescale(as.data.frame(apply(stereo1,2,s.content3, itemLevel=TRUE)))
content2 <- rescale(as.data.frame(apply(stereo2,2,s.content3, itemLevel=TRUE)))

## Study 1
content1$id <- seq_len(nrow(content1))

B <- 1500
set.seed(43215)
content.boot1 <- array(NA, dim=list(ncol(content1)-1,B)) #Variable, bootstraps
for (i in seq_len(B)){
	for (j in 1:28){ #Variables
		resample <- sample(content1$id, nrow(content1), replace=T)
		temp <- content1[resample,j]
		content.boot1[j,i] <- mean(temp,na.rm=TRUE)
	}
}

content.mat1 <- data.frame(V1=apply(content.boot1,1,mean,na.rm=TRUE), V2=apply(content.boot1,1,quantile,0.025,na.rm=TRUE), V3=apply(content.boot1,1,quantile,0.975,na.rm=TRUE), V4=c(1:28))
content.mat1$V5 <- "Domestic"
content.mat1$V5[c(1:12,13:14,17:20,21:22)] <- "Foreign" 
content.mat1$V6 <- factor(c("Stop piracy", "Stop piracy (unilateral)", "Stop piracy (multilateral)", "Enviro sanction", "Enviro sanction (unilateral)", "Enviro sanction (multilateral)", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Troops humanitarian", "Troops humanitarian (unilateral)", "Troops humanitarian (multilateral)", "Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion"), levels=c("Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops humanitarian", "Troops humanitarian (unilateral)", "Troops humanitarian (multilateral)", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Enviro sanction", "Enviro sanction (unilateral)", "Enviro sanction (multilateral)", "Stop piracy", "Stop piracy (unilateral)", "Stop piracy (multilateral)"))
content.mat1$V7 <- content.mat1$V5
content.mat1$V7[1:12] <- "Foreign (Approach)"
content.mat1$V8 <- c(rep(1:4,each=3), rep(5:12,each=2))

dev.new(height=7, width=8.8)
ggplot(data=content.mat1, aes(x=V1,y=V6)) + geom_point() + geom_segment(aes(x=V2, y=V6,xend=V3, yend=V6)) +theme_bw() + guides(color=FALSE, alpha=FALSE) + scale_color_grey() + facet_wrap(~V7, nrow=3, scales="free") +  labs(x="Stereotype content", y="") + geom_line(aes(x=V1, y=V6, group=V8, alpha=0.8)) + xlim(0,1)
#Edit group ordering in the third panel manually

##### Figure 13: Stereotype content dropping "neither" category, study 2

content2$id <- seq_len(nrow(content2))

B <- 1500
set.seed(43215)
content.boot2 <- array(NA, dim=list(ncol(content2)-1,B)) #Variable, bootstraps
for (i in seq_len(B)){
	for (j in 1:23){ #Variables
		resample <- sample(content2$id, nrow(content2), replace=T)
		temp <- content2[resample,j]
		content.boot2[j,i] <- mean(temp,na.rm=TRUE)
	}
}

content.mat2 <- data.frame(V1=apply(content.boot2,1,mean,na.rm=TRUE), V2=apply(content.boot2,1,quantile,0.025,na.rm=TRUE), V3=apply(content.boot2,1,quantile,0.975,na.rm=TRUE), V4=c(1:23))
content.mat2$V5 <- "Domestic"
content.mat2$V5[c(1:3,4:5,8:13,22:23)] <- "Foreign" 
content.mat2$V6 <- factor(c("Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)","Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion", "Anti-affirmative action", "Pro-affirmative action", "Anti-immigration", "Pro-immigration"), levels=c("Pro-affirmative action", "Anti-affirmative action", "Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-immigration", "Anti-immigration", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)"))
content.mat2$V7 <- content.mat2$V5
content.mat2$V7[1:3] <- "Foreign (Approach)"
content.mat2$V8 <- c(rep(1,3), rep(2:11,each=2))

dev.new(height=7, width=8.8)
ggplot(data=content.mat2, aes(x=V1,y=V6)) + geom_point() + geom_segment(aes(x=V2, y=V6,xend=V3, yend=V6)) +theme_bw() + guides(color=FALSE, alpha=FALSE) + scale_color_grey() + facet_wrap(~V7, nrow=3, scales="free") +  labs(x="Stereotype content", y="") + geom_line(aes(x=V1, y=V6, group=V8, alpha=0.8)) + xlim(0,1)
#Edit group ordering in the third panel manually

##### Figure 14: Stereotype content by "neither" classification rates

neither.plot <- data.frame(x=c(apply(stereo1, 2, s.contentN, itemLevel=FALSE), apply(stereo2, 2, s.contentN, itemLevel=FALSE)), y=c(content.mat1$V1, content.mat2$V1), for.lab=c(for.lab1,for.lab2))

dev.new(height=4,width=5.25)
ggplot(neither.plot, aes(x,y,color=for.lab)) + geom_point() + theme_bw() + labs(x="Pr(Neither)", y="Stereotype content (amended midpoint)") + ylim(0,1) + xlim(0,1) + guides(color=guide_legend(title="Issue type"))

#Calculate difference between foreign and domestic issues
t.test(neither.plot$x ~ neither.plot$for.lab)

#What's the one rogue issue? Isolationism
which(neither.plot$x > 0.24) 

#What about isolationism in 2018?
neither.plot$x[28+8] #Drops to 17%


####### Appendix 2.3: Political sophistication and stereotype accuracy

#### Table 6: More politically sophisticated respondents are more aware of polarization and see foreign policy issues as relatively less stereotypical

mod.1 <- lmer(value ~ 1 + age + male + white + age + educ1 + pid1 + polInterest1 + Foreign + Foreign*educ1 + (1|id) + (1|issue) + (1|year), data=types)

mod.2 <- lmer(value ~ 1 + age + male + white + age + educ1 + pid1 + polInterest1 + Foreign + Foreign*polInterest1 + (1|id) + (1|issue) + (1|year), data=types)

mod.3 <- lmer(value ~ 1 + age + male + white + age + educ1 + pid1 + polInterest1 + Polarization + Polarization*educ1 + (1|id) + (1|issue) + (1|year), data=types)

mod.4 <- lmer(value ~ 1 + age + male + white + age + educ1 + pid1 + polInterest1 + Polarization + Polarization*polInterest1 + (1|id) + (1|issue) + (1|year), data=types) 

stargazer(mod.1, mod.2, mod.3, mod.4, title="Polarization and stereotypicality", omit.stat=c("LL", "ser", "f"), style="apsr", digits=3, label="tab:stereotypicality2", covariate.labels=c("Age", "Male", "White", "Education", "Partisanship", "Political interest", "Foreign policy issue", "Foreign policy issue x Education", "Foreign policy issue x Political interest", "Polarization", "Polarization x Education", "Polarization x Political interest", "Constant"))

#### Figure 15: Political sophistication interaction plots

#Left-hand panel: Education as moderator
Bz <- fixef(mod.1)[5]
Bx <- fixef(mod.1)[8]
Bxz <- fixef(mod.1)[9]
educ.x <- seq(0,1,length.out=10)
v.mod <- vcov(mod.1)[8,8] + 2*educ.x*vcov(mod.1)[8,9] + (educ.x^2)*vcov(mod.1)[9,9] #Variance

effect.mod <- matrix(0,length(educ.x),5)
effect.mod[,1] <- Bx + Bxz*educ.x
effect.mod[,2] <- effect.mod[,1] +1.96*sqrt(v.mod)
effect.mod[,3] <- effect.mod[,1] +1.64*sqrt(v.mod)
effect.mod[,4] <- effect.mod[,1] -1.64*sqrt(v.mod)
effect.mod[,5] <- effect.mod[,1] -1.96*sqrt(v.mod)

dev.new(height=5.1, width=5)
par(oma=c(0,0,0,0), mar=c(4.1,4.1,2,1))
plot(0,0,type="n", xlim=c(0,1), ylim=c(-0.3,0), axes=FALSE, ylab="Conditional effect of foreign policy issue", xlab="Education")
polygon(x=c(educ.x, rev(educ.x)), y=c(effect.mod[,2], rev(effect.mod[,5])), col=rgb(0,0,0,0.3), border=NA)
lines(educ.x,effect.mod[,1], lwd=1.5)
abline(h=0, lty=3)
axis(2)
axis(1, at=c(0,0.5,1))
box()

#Right-hand panel: Political interest as moderator
Bz <- fixef(mod.2)[7]
Bx < fixef(mod.2)[8]
Bxz <- fixef(mod.2)[9]
int.x <- seq(0,1,length.out=10)
v.mod <- vcov(mod.1)[8,8] + 2*int.x*vcov(mod.1)[8,9] + (int.x^2)*vcov(mod.1)[9,9] #Variance

effect.mod <- matrix(0,length(int.x),5)
effect.mod[,1] <- Bx + Bxz*int.x
effect.mod[,2] <- effect.mod[,1] +1.96*sqrt(v.mod)
effect.mod[,3] <- effect.mod[,1] +1.64*sqrt(v.mod)
effect.mod[,4] <- effect.mod[,1] -1.64*sqrt(v.mod)
effect.mod[,5] <- effect.mod[,1] -1.96*sqrt(v.mod)

dev.new(height=5.1, width=5)
par(oma=c(0,0,0,0), mar=c(4.1,4.1,2,1))
plot(0,0,type="n", xlim=c(0,1), ylim=c(-0.3,0), axes=FALSE, ylab="Conditional effect of foreign policy issue", xlab="Political interest")
polygon(x=c(int.x, rev(int.x)), y=c(effect.mod[,2], rev(effect.mod[,5])), col=rgb(0,0,0,0.3), border=NA)
lines(int.x,effect.mod[,1], lwd=1.5)
abline(h=0, lty=3)
axis(2)
axis(1, at=c(0,0.5,1))
box()

##### Table 7: Ordinal cumulative link mixed-effect model

types$valueFactor <- as.factor(types$value)
mod.1o <- clmm(valueFactor ~ year + (1|id) + (1|issue), data=types, link="probit")

mod.2o <- clmm(valueFactor ~ 1 + age + male + white + age + educ1 + pid1 + polInterest1 + year + (1|id) + (1|issue), data=types, link="probit")

mod.3o <- clmm(valueFactor ~ 1 + age + male + white + age + educ1 + pid1 + polInterest1 + prefsStr + year + (1|id) + (1|issue), data=types, link="probit")

mod.4o <- clmm(valueFactor ~ 1 + age + male + white + age + educ1 + pid1 + polInterest1 + Foreign + year + (1|id) + (1|issue), data=types, link="probit")

mod.5o <- clmm(valueFactor ~ 1 + age + male + white + age + educ1 + pid1 + polInterest1 + Polarization + year + (1|id) + (1|issue), data=types, link="probit")

#Stargazer doesn't work with CLMM objects, so clunkily generate your own model output

tab.order <- c("year2018", "age", "male", "white", "educ1", "pid1", "polInterest1", "prefsStr", "Foreign", "Polarization", "0|0.5", "0.5|1")

coef.tab <- matrix(NA, nrow=length(tab.order), ncol=5)
coef.tab[match(names(coef(mod.1o)), tab.order),1] <- coef(mod.1o)
coef.tab[match(names(coef(mod.2o)), tab.order),2] <- coef(mod.2o)
coef.tab[match(names(coef(mod.3o)), tab.order),3] <- coef(mod.3o)
coef.tab[match(names(coef(mod.4o)), tab.order),4] <- coef(mod.4o)
coef.tab[match(names(coef(mod.5o)), tab.order),5] <- coef(mod.5o)

se.tab <- matrix(NA, nrow=length(tab.order), ncol=5)
se.1 <- sqrt(diag(vcov(mod.1o)))
se.tab[na.omit(match(names(se.1)[1:(length(se.1)-2)], tab.order)),1] <- se.1[1:(length(se.1)-2)]
se.2 <- sqrt(diag(vcov(mod.2o)))
se.tab[na.omit(match(names(se.2)[1:(length(se.2)-2)], tab.order)),2] <- se.2[1:(length(se.2)-2)]
se.3 <- sqrt(diag(vcov(mod.3o)))
se.tab[na.omit(match(names(se.3)[1:(length(se.3)-2)], tab.order)),3] <- se.3[1:(length(se.3)-2)]
se.4 <- sqrt(diag(vcov(mod.4o)))
se.tab[na.omit(match(names(se.4)[1:(length(se.4)-2)], tab.order)),4] <- se.4[1:(length(se.4)-2)]
se.5 <- sqrt(diag(vcov(mod.5o)))
se.tab[na.omit(match(names(se.5)[1:(length(se.5)-2)], tab.order)),5] <- se.5[1:(length(se.5)-2)]

z.tab <- coef.tab/se.tab
p.tab <- 2*pnorm(abs(z.tab), lower.tail=FALSE)
star.tab <- ifelse(p.tab<0.001, "***", ifelse(p.tab<0.05, "**", ifelse(p.tab<0.1,"*","")))
coef.tab <- round(coef.tab, digits=3)
se.tab <- round(se.tab, digits=3)

out.tab <- matrix(NA, nrow=26, ncol=5)
out.tab[seq(1,23,by=2),] <- paste(coef.tab,star.tab, sep="")
out.tab[seq(2,24,by=2),] <- paste("(",se.tab,")", sep="")
out.tab[which(out.tab=="(NA)" | out.tab=="NANA")] <- NA
out.tab[25,] <- c(mod.1o$dims$n, mod.2o$dims$n, mod.3o$dims$n, mod.4o$dims$n, mod.5o$dims$n)
out.tab[26,] <- round(c(BIC(mod.1o), BIC(mod.2o), BIC(mod.3o), BIC(mod.4o), BIC(mod.5o)),digits=2)
rownames(out.tab) <- rep("",26)
rownames(out.tab)[seq(1,23,by=2)] <- tab.order
rownames(out.tab)[25:26] <- c("N","BIC")
out.tab

stargazer(out.tab, title="Ordinal cumulative link mixed-effect model", style="apsr", digits=3, label="tab:ordinal")

####### Appendix 2.4 Order effects

### Figure 16a: Marginal effects of order on intensity

#Obtain question sequence ordering for the policy preference series (series 1) and the stereotype series (series 2)
seq1 <- c("Q4", "Q5", "Q6", "Q7", "Q8", "Q9", "Q10", "Q11", "Q12", "Q13", "Q14", "Q15") #The question ID sequence for Series 1
seq2 <- c("Q17", "Q18", "Q19", "Q20", "Q21", "Q22", "Q23", "Q24", "Q25", "Q26", "Q27", "Q28") #The question ID sequence for Series 2

#Convert ordering to matrix
order1.mat <- strsplit(dat1$DO.BL.QuestionSeries1, split="\\|")
order2.mat <- strsplit(dat1$DO.BL.QuestionSeries2, split="\\|")

#Examples of how it works
#match(seq1, order1.mat[[1]]) - 1 #Q4 appears 5th, Q5 appears 12th, etc. (we remove a constant, as the first question is just an intro)
#match(seq2, order2.mat[[1]]) - 1 #Q17 appears 12th, Q18 appears 1st, etc. (we remove a constant, as the first question is just an intro)

reorder1.list <- reorder2.list <- list(NA)
orderdiff <- matrix(nrow=nrow(dat1), ncol=length(seq1))
for (i in 1:length(order1.mat)){
	reorder1.list[[i]] <- match(seq1, order1.mat[[i]]) - 1
	reorder2.list[[i]] <- match(seq2, order2.mat[[i]]) - 1 
	orderdiff[i,] <- reorder2.list[[i]] + (12-reorder1.list[[i]])
}
orderdiff2 <- rescale(orderdiff[,c(rep(1:4,each=3),rep(5:12,each=2))])

coefs <- matrix(NA, nrow=28, ncol=2)
for (i in 1:28){
	mod <- lm(intensity1[,i] ~ orderdiff2[,i] + policy1[,i] + dat1$pid1 + dat1$educ1 + dat1$male + dat1$white + dat1$age + dat1$polInterest1 + dat1$income)
	coefs[i,1] <- coef(mod)[2]
	coefs[i,2] <- sqrt(vcov(mod)[2,2])
}

dev.new(height=5,width=5.2)
par(oma=c(0,0,0,0), mar=c(2,4.1,2,2))
plot(0,0,type="n", xlim=c(1,28), ylim=c(-1,1), axes=FALSE, ylab="Effect on stereotype intensity", xlab="")
points(1:28,coefs[,1], pch=16)
for (i in 1:28){
	lines(x=c(i,i), y=c(coefs[i,1]-1.96*coefs[i,2], coefs[i,1]+1.96*coefs[i,2]))
}
abline(h=0)
axis(2)
box()

### Figure 16b: Conditional effects of order on intensity

policy1.s <- abs(policy1[,1:28]-4)/3

coefs2 <- matrix(NA, nrow=28, ncol=2)
for (i in 1:28){
	mod <- lm(intensity1[,i] ~ orderdiff2[,i]*policy1.s[,i] + dat1$pid1 + dat1$educ1 + dat1$male + dat1$white + dat1$age + dat1$polInterest1 + dat1$income)
	coefs2[i,1] <- coef(mod)[11]
	coefs2[i,2] <- sqrt(vcov(mod)[11,11])
}
dev.new(height=5,width=5.2)
par(oma=c(0,0,0,0), mar=c(2,4.1,2,2))
plot(0,0,type="n", xlim=c(1,28), ylim=c(-1,1), axes=FALSE, ylab="Effect on stereotype intensity", xlab="")
points(1:28,coefs2[,1], pch=16)
for (i in 1:28){
	lines(x=c(i,i), y=c(coefs2[i,1]-1.96*coefs2[i,2], coefs2[i,1]+1.96*coefs2[i,2]))
}
abline(h=0)
axis(2)
box()

### Figure 16c: Marginal effects of order on stereotype content

coefs3 <- matrix(NA, nrow=28, ncol=2)
for (i in 1:28){
	mod <- lm(content1[,i] ~ orderdiff2[,i] + policy1[,i] + dat1$pid1 + dat1$educ1 + dat1$male + dat1$white + dat1$age + dat1$polInterest1 + dat1$income)
	coefs3[i,1] <- coef(mod)[2]
	coefs3[i,2] <- sqrt(vcov(mod)[2,2])
}

dev.new(height=5,width=5.2)
par(oma=c(0,0,0,0), mar=c(2,4.1,2,2))
plot(0,0,type="n", xlim=c(1,28), ylim=c(-1,1), axes=FALSE, ylab="Effect on stereotype content", xlab="")
points(1:28,coefs3[,1], pch=16)
for (i in 1:28){
	lines(x=c(i,i), y=c(coefs3[i,1]-1.96*coefs3[i,2], coefs3[i,1]+1.96*coefs3[i,2]))
}
abline(h=0)
axis(2)
box()

### Figure 16d: Conditional effects of order on stereotype content

coefs4 <- matrix(NA, nrow=28, ncol=2)
for (i in 1:28){
	mod <- lm(content1[,i] ~ orderdiff2[,i]*policy1.s[,i] + dat1$pid1 + dat1$educ1 + dat1$male + dat1$white + dat1$age + dat1$polInterest1 + dat1$income)
	coefs4[i,1] <- coef(mod)[11]
	coefs4[i,2] <- sqrt(vcov(mod)[11,11])
}

dev.new(height=5,width=5.2)
par(oma=c(0,0,0,0), mar=c(2,4.1,2,2))
plot(0,0,type="n", xlim=c(1,28), ylim=c(-1,1), axes=FALSE, ylab="Effect on stereotype content", xlab="")
points(1:28,coefs4[,1], pch=16)
for (i in 1:28){
	lines(x=c(i,i), y=c(coefs4[i,1]-1.96*coefs4[i,2], coefs4[i,1]+1.96*coefs4[i,2]))
}
abline(h=0)
axis(2)
box()

###### Appendix 2.5: Cohort effects

### Table 8: GAMs
 
#Generate average stereotype intensity measures 
dat1$hawkStereo <- apply(intensity1[,17:22],1,mean,na.rm=TRUE) #Hawkishness/interventionism
dat1$fissStereo <- apply(intensity1[,c(13:14,17:22)],1,mean,na.rm=TRUE) #Foreign policy issues
dat1$fappStereo <- apply(intensity1[,1:12],1,mean,na.rm=TRUE) #Foreign policy approaches
dat1$fpStereo <- apply(intensity1[,c(1:12,13:14,17:22)],1,mean,na.rm=TRUE) #All foreign issues
dat1$domStereo <- apply(intensity1[,c(15:16,23:28)],1,mean,na.rm=TRUE) 
dat1$fdStereo <- dat1$fpStereo-dat1$domStereo #Foreign - domestic

mod.1 <- gam(hawkStereo ~ s(age, bs="cr") + male + income + white + educ + pid + polInterest, data=dat1)
mod.2 <- gam(fissStereo ~ s(age, bs="cr") + male + income + white + educ + pid + polInterest, data=dat1)
mod.3 <- gam(fappStereo ~ s(age, bs="cr") + male + income + white + educ + pid + polInterest, data=dat1)
mod.4 <- gam(fpStereo ~ s(age, bs="cr") + male + income + white + educ + pid + polInterest, data=dat1)
mod.5 <- gam(fdStereo ~ s(age, bs="cr") + male + income + white + educ + pid + polInterest, data=dat1)

stargazer(mod.1, mod.2, mod.3, mod.4, mod.5, title="Generalized additive models find little evidence of cohort effects", omit.stat=c("LL", "ser", "f"), style="apsr", digits=3, label="gam", covariate.labels=c("Male", "Income", "White", "Education", "Partisanship", "Pol. Interest"))
#Manually add statistics for spline on age:

#EDF, Chisq: 
#4.586, 1.944*
#4.869, 1.688
#5.111, 4.563***
#5.47, 3.331***
#2.261, 4.333***

#### Figure 17: Cohort effects

#Panel a: hawkishness
p <- plot(mod.1, select=1, rug=TRUE, se=TRUE, ylab="Average stereotype intensity", residual=TRUE, xlab="Age")

dev.new(height=5,width=5.2)
par(mar=c(4,5,0.5,0.5))
plot(0,0, xlim=c(18,95), ylim=c(-0.7,0.7), type="n", xlab="Age", ylab="Average stereotype intensity:\n hawkish/interventionist policies")
points(p[[1]]$raw, p[[1]]$p.resid, pch=16, cex=0.5, col=rgb(0,0,0,0.25))
lines(p[[1]]$x, p[[1]]$fit)
polygon(x=c(p[[1]]$x,rev(p[[1]]$x)),y=c(p[[1]]$fit+p[[1]]$se,rev(p[[1]]$fit-p[[1]]$se)), col=rgb(0,0,0,0.25), border=rgb(1,1,1,0.25))

#Panel b: Foreign policies
p <- plot(mod.2, select=1, rug=TRUE, se=TRUE, ylab="Average stereotype intensity", residual=TRUE, xlab="Age")

dev.new(height=5,width=5.2)
par(mar=c(4,5,0.5,0.5))
plot(0,0, xlim=c(18,95), ylim=c(-0.7,0.7), type="n", xlab="Age", ylab="Average stereotype intensity:\n foreign policy issues")
points(p[[1]]$raw, p[[1]]$p.resid, pch=16, cex=0.5, col=rgb(0,0,0,0.25))
lines(p[[1]]$x, p[[1]]$fit)
polygon(x=c(p[[1]]$x,rev(p[[1]]$x)),y=c(p[[1]]$fit+p[[1]]$se,rev(p[[1]]$fit-p[[1]]$se)), col=rgb(0,0,0,0.25), border=rgb(1,1,1,0.25))

#Panel c: Foreign policy approaches
p <- plot(mod.3, select=1, rug=TRUE, se=TRUE, ylab="Average stereotype intensity", residual=TRUE, xlab="Age")

dev.new(height=5,width=5.2)
par(mar=c(4,5,0.5,0.5))
plot(0,0, xlim=c(18,95), ylim=c(-0.7,0.7), type="n", xlab="Age", ylab="Average stereotype intensity:\n foreign policy approaches")
points(p[[1]]$raw, p[[1]]$p.resid, pch=16, cex=0.5, col=rgb(0,0,0,0.25))
lines(p[[1]]$x, p[[1]]$fit)
polygon(x=c(p[[1]]$x,rev(p[[1]]$x)),y=c(p[[1]]$fit+p[[1]]$se,rev(p[[1]]$fit-p[[1]]$se)), col=rgb(0,0,0,0.25), border=rgb(1,1,1,0.25))

#Panel d: All foreign policy issues
p <- plot(mod.4, select=1, rug=TRUE, se=TRUE, ylab="Average stereotype intensity", residual=TRUE, xlab="Age")

dev.new(height=5,width=5.2)
par(mar=c(4,5,0.5,0.5))
plot(0,0, xlim=c(18,95), ylim=c(-0.7,0.7), type="n", xlab="Age", ylab="Average stereotype intensity:\n foreign policy issues and approaches")
points(p[[1]]$raw, p[[1]]$p.resid, pch=16, cex=0.5, col=rgb(0,0,0,0.25))
lines(p[[1]]$x, p[[1]]$fit)
polygon(x=c(p[[1]]$x,rev(p[[1]]$x)),y=c(p[[1]]$fit+p[[1]]$se,rev(p[[1]]$fit-p[[1]]$se)), col=rgb(0,0,0,0.25), border=rgb(1,1,1,0.25))

#Panel e: Foreign - domestic
p <- plot(mod.5, select=1, rug=TRUE, se=TRUE, ylab="Average stereotype intensity", residual=TRUE, xlab="Age")

dev.new(height=5,width=5.2)
par(mar=c(4,5,0.5,0.5))
plot(0,0, xlim=c(18,95), ylim=c(-0.7,0.7), type="n", xlab="Age", ylab="Average stereotype intensity:\n foreign policy - domestic policy")
points(p[[1]]$raw, p[[1]]$p.resid, pch=16, cex=0.5, col=rgb(0,0,0,0.25))
lines(p[[1]]$x, p[[1]]$fit)
polygon(x=c(p[[1]]$x,rev(p[[1]]$x)),y=c(p[[1]]$fit+p[[1]]$se,rev(p[[1]]$fit-p[[1]]$se)), col=rgb(0,0,0,0.25), border=rgb(1,1,1,0.25))

###### Appendix 2.6: Weighted analysis

### Figure 18: Changes in second-order beliefs from 2014-18 reflect changes in first-order preferences (weighted)
stereo.dat$weight[which(is.na(stereo.dat$weight))] <- 0
temp.18 <- which(stereo.dat$year==2018)
temp.14 <- which(stereo.dat$year==2014)

dif.content <- rep(NA,19)
for (j in 1:19){ #Variables
	dif.content[j] <- weighted.mean(stereo.dat[temp.18,j],stereo.dat$weight[temp.18],na.rm=TRUE) - mean(stereo.dat[temp.14,j],na.rm=TRUE)
	}

dif.policy <- rep(NA,19)
policy.dat$weight[which(is.na(policy.dat$weight))] <- 0
for (j in 1:19){ #Variables
	temp.18.d <- which(policy.dat$year==2018 & policy.dat$pid3==1)
	temp.14.d <- which(policy.dat$year==2014 & policy.dat$pid3==1)
	temp.18.r <- which(policy.dat$year==2018 & policy.dat$pid3==3)
	temp.14.r <- which(policy.dat$year==2014 & policy.dat$pid3==3)
	dif.policy[j] <- (weighted.mean(policy.dat[temp.18.r,j],w=policy.dat$weight[temp.18.r], na.rm=TRUE) - mean(policy.dat[temp.14.r,j],na.rm=TRUE)) - (weighted.mean(policy.dat[temp.18.d,j],w=policy.dat$weight[temp.18.d], na.rm=TRUE) - mean(policy.dat[temp.14.d,j],na.rm=TRUE))
}

dif.plot2 <- data.frame(stereo=dif.content, policy=dif.policy, V6=factor(c("Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)", "Protectionist", "Free trade", "Anti-environment", "Pro-environment", "Isolationist", "Interventionist", "Anti-arms control", "Pro-arms control", "Decrease military spending", "Increase military spending", "Anti-gun control", "Pro-gun control", "Anti-tax hike", "Pro-tax hike", "Anti-abortion", "Pro-abortion"), levels=c("Pro-abortion", "Anti-abortion", "Pro-tax hike", "Anti-tax hike", "Pro-environment", "Anti-environment", "Pro-gun control", "Anti-gun control", "Pro-arms control", "Anti-arms control", "Interventionist", "Isolationist", "Free trade", "Protectionist", "Increase military spending", "Decrease military spending", "Troops oil", "Troops oil (unilateral)", "Troops oil (multilateral)"))) 

mod.r2 <- lm(dif.content ~ dif.policy)
summary(mod.r2)$r.squared #r2= 0.72	
cor(dif.content, dif.policy, use="complete.obs") #r=0.85

dev.new(height=7, width=7)
ggplot(dif.plot2, aes(x=policy, y=stereo)) + geom_point() + theme_bw() + guides(alpha=FALSE) + labs(x="Difference-in-difference in policy preferences: Republicans vs Democrats, 2014-2018", y="Change in stereotype content, 2014-18") + geom_vline(xintercept=0, lty=3) + geom_hline(yintercept=0, lty=3) + geom_smooth() + geom_label_repel(aes(policy,stereo,label=V6), data=dif.plot2, force=1, show.legend=FALSE, size=3) 
